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Abstract 

Two-dimensional correlation spectroscopy (2DCS) based on the nonlinear optical response of 
excitons to sequences of ultrafast pulses, has the potential to provide some unique insights into 
carrier dynamics in semiconductors. The most prominent feature of 2DCS, cross peaks, can best 
be understood using a sum-over-states picture involving the many-body eigenstates. However, the 
optical response of semiconductors is usually calculated by solving truncated equations of motion 
for dynamical variables, which result in a quasiparticle picture. In this work we derive Green's 
function expressions for the four wave mixing signals generated in various phase-matching directions 
and use them to establish the connection between the two pictures. The formal connection with 
Frenkel excitons (hard-core bosons) and vibrational excitons (soft-core bosons) is pointed out. 

PACS numbers: 78.47.+p,71.35.-y 
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I. INTRODUCTION 



Exciton models are widely used to describe the linear and nonlinear optical properties of 
many types of systems, including semiconductor nanostructures (quantum wells, dots and 
wires), molecular aggregates and crystals,-^ as well as vibrations in proteins.-^ In semi- 
conductors, nonlinear optical experiments reveal a wealth of interesting phenomena.- 1 ^"^ 
For instance, such experiments provide information about many-exciton states such as biex- 
citons, their interactions, relaxation and dissociation.— »i2iiiLM 

The introduction of multidimensional techniques had revolutionized NMR in the 
seventies^ and established it as a powerful tool for studying complex systems and iden- 
tifying specific structural and dynamical correlations.— In such experiments the system is 
subjected to a sequence of well separated pulses. Correlation plots of the signals vs. two 
(or more) time delay periods then provide multidimensional spectroscopic windows into the 
system. The correlated dynamics of spins carefully prepared by the pulse sequence is very 
sensitive to their interactions. Analysis of these correlation plots then provides a powerful 
probe for molecular geometries and dynamical correlations. These techniques were recently 
extended to the infrared and the visible regime and were shown to be very useful for Frenkel 
excitons in molecular systems.- 1 ^ 1 ^ 1 ^ There are some differences between the optical and 
the NMR techniques. NMR uses strong saturating fields whereas optical pulses are most 
effective in the weak field regime. NMR signals are essentially isotropic in space whereas 
coherent optical signals are generated in well defined (phase-matching) directions. These 



differences were explored in detail in Refs. I20l j2l|j22l . Nevertheless the NMR and optical 



techniques are conceptually similar and many ideas of pulse sequences developed in NMR 
may be adopted in the optical regime, where the millisecond NMR time-scale is pushed to 
the femtosecond regime. The same ideas may be extended to study interband and inter- 
suband excitations in semiconductors.— i2ii2£i2&2Z Multidimensional analysis of the nonlinear 
optical response of semiconductors to sequences of femtosecond pulses could provide a novel 
probe for many-body interactions. In a recent work 23 on semiconductor Quantum Wells, 
2D correlation spectra from three 3rd order optical techniques have been calculated. The 
unique character of 2D spectroscopy allowed to easily recognize and classify features due to 
different types of biexcitons. Such features are sometimes difficult to separate in the usual 
one-dimensional mode of displaying non-linear spectra, due to the strong line broadening 
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and the highly congested exciton spectra. 

Two types of approaches have been traditionally used towards modeling the nonlinear op- 
tical response of excitonic systems. The first is based on the many-body eigenstates obtained 
by exact diagonalization of the Hamiltonian.— Sum-over-states (SOS) expressions can then 
be derived for the nonlinear response functions and optical signals. This method is practical 
in many applications to electronic and vibrational Frenkel excitons in molecules^^i and 
allows clear identification and classification of possible single- and multi-photon resonances. 
Calculating the eigenstates is a serious computational bottleneck in extended structures. 
For an iV site tight-binding Frenkel-exciton model the number of single and two-exciton 
states scales as ~ N and ~ iV 2 respectively. For Wannier excitons in semiconductors these 
scalings are ~ iV 2 and ~ iV 4 , making the simulations prohibitively expensive. This is why 
the approach is not widely used for electron-hole excitations in semiconductors. Instead, 
one adopts a second strategy, which describes the response in terms of quasiparticles (QP), 
and the many-particle eigenstates are never calculated.- 1 ^'^^ 4 -'^ Calculations are per- 
formed by solving equations of motion for microscopic coherences, which are coupled to 
other dynamical variables. Even for a simple system such as a single semiconductor quan- 
tum well, solving the equations numerically to create a 2D map of a nonlinear response 
function is computationally expensive,— since these equations must be solved repeatedly for 
different pulse delays. Only after obtaining the optical signal on a 2D time grid, a Fourier 
transform can be performed to get the 2DCS. Apart from direct, numerical solutions of 
equations of motion 3 ^ 1 ^ there exist other theoretical approaches to exciton correlation ef- 
fects, such as memory kernel representation 3 ^ 3 ^ or Coupled Cluster Expansion for doped 
semi conduct ors .— »i 

In this paper we derive closed expressions for 2DCS of semiconductors by solving the 
Nonlinear Exciton Equations (NEE) 3,42 for the third order response. Both time-ordered 
and non-ordered forms of the response function which represent time and frequency domain 
techniques, respectively, are derived. Our QP expressions for the response are given in terms 
of the single exciton Green's function and the exciton scattering matrix. The SOS response 
functions, in contrast, are expressed in terms of many-exciton eigenstates. Even though the 
response functions calculated using both techniques must be identical, the relation between 
the two pictures is not obvious. The expressions look very different and it is not possible to 
see their equivalence by a simple inspection. The SOS expressions contain large terms, which 
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grow with system size and have opposite signs, thus they almost cancel. This complicates 
their numerical implementation. In contrast these cancellations are built-in from the outset 
in the QP approach, which uses a harmonic reference system. The nonlinearities are then 
attributed to exciton-exciton scattering which is absent in the harmonic reference system. 
The second goal of this paper is to show precisely how the two pictures of many-body 
correlations are connected. We write down the SOS expressions using the Keldysh loop 
and then derive the QP expressions directly from the SOS ones. This provides a time- 
domain interpretation for the interference effects. The SOS and the QP expressions provide 
complementary views into the origin of features seen in 2D spectrograms. 

In Sec. |TT] we present the SOS expressions for the third order response obtained from 
time-dependent perturbation theory. Their QP counterparts are derived in Sec. IHIi We 
use the method developed in Refs. [3] 34 



to transform the Hamiltonian to a form typical 
for interacting oscillators. The starting many-electron Hamiltonian can be written in an 
ab-initio,— tight-binding 44 or a k • p basis. One of the key results of this paper, i.e., the 
equivalence of the SOS and QP pictures is proven in Sec. [TV] In Sec. [V] we derive closed 
expressions for 2D correlation signals. The QP approach provides a unified description 
for electron-hole excitations in semiconductors as well as to Frenkel excitons in molecular 
aggregates (Paulions) and anharmonic vibrations (bosons), which are described by the same 
general Hamiltonian. QP formulae for nonlinear response have been derived previously 
along similar lines for Frenkel excitons. This connection is shown in Appendix [Fl In the 
last Section ( jVIl) we discuss the results. 



II. SUM-OVER-STATES EXPRESSIONS FOR THE TIME-ORDERED NONLIN- 
EAR RESPONSE 

We consider a 4 wave-mixing experiment performed with three femtosecond laser pulses 
(Fig. [T|). The optical electric field is: 

3 

E (r, t) = J2 E, (r, t) = E + (r, t) + E~(r, t), (1) 
i=i 

3 

E + (r,t) = ^S+{t - r j )e- iu t t e ik * T , (2) 
i=i 
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3 

E~(r,t) = J2 S i^ ~ 7j)e^*e-^ r . (3) 

3=1 

The j-th pulse is centered at Tj, has an envelope Sj(t — Tj), carrier frequency Uj, and 
wavevector k,. E + (E~~) denotes the positive (negative) frequency part of the field, and 
£j = • The induced polarization in the system is recorded as a function of time- 

delays between pulses. 

Assuming the dipole interaction with the optical field Hj = fi ■ E(r,r), where ft is the 
dipole operator, the third-order contribution to the system's polarization can be written as 



P(r,r i )= /// dT 3 dT 2 dT 1 S^ sos \n,T 3 ,T 2 ,r 1 )E(r,r 3 )E(r,r 2 )E(r,r 1 ), (4) 




where the response function S^ sos \ which connects the induced polarization with the laser 
field envelopes, is given by (throughout this paper we set K = 1): 

5( 5 ° s )(r 4 ,r 3 ,r 2 ,r 1 ) = t 3 [0(t 43 )0(t 32 )0(t 21 ) (fx(r 4 ) p,(r 3 ) p,(r 2 ) p,( Tl )) (5) 
- 6(t 43 )6(t A2 )6(t 21 ) (A(t 3 )A(t4)A(t2)A(ti)) 
+ d(r 42 )9(r 23 )9(T 41 ) (fi(r 3 ) fi(r 2 ) fi(r 4 ) 
-e(r 41 )6(r 12 )6(r 23 ) (frfo) fifa) fi(n) ftn))] . 

We shall use double-sided Feynman diagrams to represent the time ordering of various 
interactions.— The four terms in Eq. (jHJ) are represented by diagrams a, b, c, d shown on 
Figure El These diagrams should be read starting at the bottom left and proceeding along 
the loop, clockwise, as indicated by the arrows. The Tj variables are ordered on the Keldysh- 
Schwinger loop, but not necessarily in real (physical) time. Tj in diagrams (a) and (d) are 
also ordered in real time. This is not the case for diagrams (b) and (c): in (b) r 3 can come 
either before or after T\ and t 2 , whereas in (c) T\ can come either before or after r 3 and r 2 . 

If the eigenstates \a) and eigenvalues e a of the system are known, Eq. (JHJ) may be expanded 
in terms of the corresponding matrix elements: 

(^(r 4 )^(r 3 )/x(r 2 )/x(r 1 )) (6) 

— V~" II II II II P~* \(£a 3 -e g )T4, + (ea 2 -ea 3 )T3+(Sa 1 -Sa 2 )T2+(Sg-Sa 1 )Tl\ 

/ , t^ga-AH 1 0.30,2 t L a,2a\ H'aig^ L J • 

(11,(12,(13 

So far we considered a general multilevel system. We next turn to the response of excitons, 
where the energy levels form manifolds, classified by the number of excitons: the ground 



state (g), single exciton (e), two-exciton (/) (or biexciton), etc. (Fig. [3]). We shall assume 
that the dipole operator can only create and annihilate a single exciton at a time. Only 
the single and the two-exciton states then contribute to the third order signals. We further 
partition the dipole operator as A — fi + + A j where A + is the positive frequency part 
which induces upward g to e and e to / transitions, while its Hermitian conjugate A (the 
negative frequency part) induces the opposite transitions. We thus write 

£„>£„/ 

a = 2J /w w) <yi • 

Invoking the rotating- wave approximation (RWA), we neglect all terms where at least one 
of the transitions is not in resonance with one of the incident carrier frequencies. The 
system-field interaction term then becomes 

Hj (t) = -fi + E + {r, r) - iTE~{r, r) 

Each correlation function in Eq. (jSJ) will split into 2 4 = 16 terms upon substituting 
A = p, + + ft ~. Assuming that the system is initially in the ground state, only two of these 
contributions are non-zero 

(AAAA) = (A~A + A~A + ) + (A~A~A + A + ) • (7) 

Substitution of Eq. (JTj) into Eq. ([3]) gives 

S^in^r^n) = * 3 [e(r 43 )9(T 32 )0(r 21 ) (^(^^(r^+ir^fi+in)) (al) (8) 
+ e(r 43 )6(r 32 )6(r 21 ) (/*" (r 4 ) A + (r 3 ) A~ (r 2 ) A + fa)) (a2) 
- 0fa3)0fa 2 )0(r 21 ) (A-fa)A~fa)A + fa)A + fa)) (b) 
+0fa 2 )0fa 3 )0(T 41 ) (A-fa)A + fa)A~fa)A + fa))] (c) 
+ c.c. 

The four terms represented by the diagrams in Fig. H] were obtained by taking A( r 4) = 
A~( r 4) f° r the last interaction, Afa) — A + ( r 4) gives the complex conjugates. Hereafter 
left /right direction of the arrows corresponds to A~/A + i n Eq.©. Note that time- reversal 
symmetry implies (A~fa)A~fa)A~fa)A~fa))* = (A + fa)A + fa)A + fa)A + 04))- If the 
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pulse envelopes are much shorter than their delays, the system is forced to interact se- 
quentially first with pulse fei, then k 2 and finally k 3 . This means that in the integral of 
Eq. (jlj) one must replace E(r, Tj) with one of the Sj, depending on the time-ordering of the 
integration variables in real (physical) time. We note that the first and the second terms 
in Eq. (|S]) impose a full time ordering of the integration variables while the third and the 
fourth terms do not. Term (b) is only partially time ordered. Depending on the position of 
r 3 relative to the T\ < r 2 < r 4 sequence, the diagram can be separated into three fully time 
ordered terms: r 3 < t 4 , t\ < r 3 < r 2 or r 2 < r 3 < r 4 . Formally we do that by separating the 
product of step functions as follows: 

^43)^42)^21) = e(u 2 )6(r 21 )d(r 13 ) + #(t 42 )£(t 23 )#(t 31 ) + 9{r i3 )e{r 32 )9{r 2l ). 

Using this relation, diagram (b) of Fig. H]is split into (b3), (b2) and (bl) as shown in the 
first line of Fig. [51 The interactions on the l.h.s. of this diagrammatic equation are ordered 
on the loop. On the other hand, the arrows in the open, double-sided diagrams on the r.h.s. 
are ordered in real (physical) time. All diagrams on the r.h.s. are obtained from (b) by 
moving the arrows while preserving their order along the loop (but not in physical time!). 
Similarly we write for term (c) 

fl(r 23 Mr 42 )#(T 41 ) = 9(t, 2 )0(t 21 )6(t 13 ) + 6(t 42 )6(t 23 )0(t 31 ) + 6(t 41 )6(t 12 )6(t 23 ) 

and the diagram is split into (c2), (c3), (cl). can now be recast in the fully time-ordered 
form 



= * 3 [9(r i3 )9(r 32 )9(T 21 ) 


(A (r 4 )A (r 3 )fi + (r 2 


)A + (n)> 


(al) (9) 


+ 0{t, 3 )6{t 32 )6{t 21 ){v 


_ (r 4 )A + (r 3 )A" 


\r 2 )^ 


(ri)> 


(a2) 


-6(T43)6(T 32 )6(T 21 )(fl- 


"(r 3 )A"(T-4)A 4 


> 2 )A + 




(bl) 


-0(T42)e(r 23 )e(r 31 )(fi- 


"(r 3 )A~(r 4 )A 4 


> 2 )A + 


Xn)) 


(b2) 


-0(t42)0(t 2 i)0(t 13 )(A~ 


"( 7 "3)A _ (^4)A 4 


-(r 2 )A + 


"(n)> 


(b3) 


+ 0(t4i)0(t 12 )0(t 23 )(A" 


"(r 3 )A + (r 2 )A" 


-(T 4 )A + (n)> 


(cl) 


+ e(T 42 )6(T 23 )6(T 31 )(fi- 


_ ( t 3)A + (t-2)A" 


-(r 4 )A + 


"(n)> 


(c2) 


+0(t42)0(t 2 i)0(t 13 )(A" 


(r 3 )A + (r 2 )A" 


(r 4 )A + 


m 


(c3) 



+ c.c. 
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The labels on the right correspond to the various diagrams shown in Figs. (j4]) and (jSj). 

Once split into fully time-ordered contributions, it is convenient to change the integra- 
tion variables in Eq. (J3J) from t±,t 3 ,t 2 ,Ti that label the actual interaction times with the 
fields, to the three delays t 3 ,t 2 ,t± between successive interactions. Note that the correla- 
tion functions are invariant to time translation (/x(r 4 — t)(X{t 3 — r)/i(r 2 — r)/i(r 1 — r)) = 
(/i(T 4 ) /x(r 3 ) /x(r 2 ) (i{t\)) . Eq. (@J thus assumes the form 

P{r,n) = dt 3 dt 2 dt l S^{t 3 ,t 2 ,t 1 )E 3 (r,T 4 -t 3 )E 2 (r,r A -t 2 -t 3 )E l {r,r A -t 1 -i 2 -t 3 ), 

(10) 

where ti = r 2 — tl, t 2 = r 3 — t 2 , t 3 = r 4 — r 3 . In the impulsive limit, where all pulses are 
shorter than all system's response time scales, we can substitute Eqs. ([I])-© in Eq. (flQj) and 
eliminate the time integrations. This gives 

P(r, r 4 ) = S {sos \t 3 , t 2 , ti)^ 3 A3 (r 4 - t 3 )£ 2 A2 (r 4 -t 3 - t 2 )E^(r 4 -t 3 -t 2 - h)x (11) 

e j(Aifci+A2fc2+A3fc3)r e -'t(Aiwi+A2W2+A3W3)r4 e i(Aiwi+A2W2+A3^3)t3 e j(Aiaji+A2W2)t2giAit<Jiti 

The polarization is created along 8 possible directions k s = X\k\ + X 2 k 2 + A 3 A; 3 with = ±1 

4 

P ( r , U ) = J2 P ( k s, w a ) e ikar - iuJaTi + c.c. (12) 

s=l 

where 

P (k s , u a ) = Si sos 1 (t 3 , t 2 , h) E^E^E^. 

k\ + k 2 + k 3 vanishes for the assumed dipole selection rules in our model. Since 
P {— k s ,— uj s ) = P*(k s ,uj s ), we are left with three independent combinations kj 
k + k 2 + k 3 , k n = +ki - k 2 + fc 3 , fc/// = +ki + k 2 - fc 3 : 

5 (5 ° s) (^,t 2 ,ti) = Sf OS) (^,t 2 ,ti) + 5}f s) (t 3 ,t 2 ,t!) + Sgf^ta^ti). 



We can classify the diagrams in Fig. according to the directions of the arrows: arrow 
pointing to the right (left) represents +k (—k), arrows are read from the bottom up on 
either side. We obtain for ki (Fig. [6]) 

Sf °^ = i 3 6(h)6(t 2 )6(t 3 ) [<A"(0)/i + (ti + * 2 )A~(*i + t2 + t 3 )A + (ti)> (c3) (13) 

+ </i-(o)/i + (ti)/i-(ti + 1 2 + t 3 )A + (ti + 1 2 )> (ci) 
- (/r (o)/r(ti + 1 2 + t 3 )A + (^i + t 2 ) / i + (t 1 )>] . (b3) 



s 



For the kjj technique we similarly have (Fig. [7]): 

s¥° s) = PefaMhMh) + h + t 3 )A + (*i + i 2 )A^i)A + (o)) (a2) (14) 

+ (A~(*i)A + (*i + t 2 )A~(*i + * 2 + t 3 )A + (o)> (c2) 

- (A~(*i)A~(*i + * 2 + t 3 )A + (ti + t 2 )A + (o)>] (b2) 

Finally fc/jj is given by (Fig. [8]): 

5^ os) = i*6{h)6{t 2 )9{h) + t 2 + * 3 )A"(*i + t 2 )A + (ti)A + (0)> (al) (15) 

- (A"(*i + h)ti~{h + t 2 + t 3 )A + (ti)A + (o)>] (bi) 



Each term is labelled according to Eq. Qj. Eqs. (I13til5j) can be used to express the third 
order SOS response in terms of transition dipoles, system frequencies and dephasing rates 
(see App. [E]and Sec. |V|) . 

In the next section we employ the EOM approach to derive the alternative QP expres- 
sions for these signals. These will then be connected with the current SOS expressions in 
Section HVl 

III. QUASIPARTICLE EXPRESSIONS FOR WANNIER EXCITONS IN SEMI- 
CONDUCTORS 

Interband transitions in semiconductors may be described by the two-band many-electron 
Hamiltonian:— ^ 

H T = H + Hc + Hj , (16) 

with the single-particle part 

1A / j t mi,ni L 'mi t 'ni ' / 4 h m 2 ,n 2 u m 2 u ™2 ) 
m\,ri\ m 2 ,n 2 

where c* create electrons and create holes. The Coulomb interaction is: 

He = ^ Kninifci/i C mi C ni C fei c 2l + 2 ^m 2 n 2 k 2 l 2 dln 2 dn 2 dk 2 di 2 

mi,ni,ki,h m 2 ,n 2 ,k 2 ,l 2 

^ 4 ^min 2 l 1 k 2 C m 1 dl l2 dk 2 Cl 1 , 
m\,n 2 ,k 2 ,l\ 
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while 

Hi = — ^ ] (£■ (t) ^ mim , 2 <^ mi d^ m2 + E (t) /i. mim2 <i m2 c mi ) , 

is the dipole interaction with light, and the optical electric field E will be treated as a scalar 
for simplicity. Ht can describe both bulk and low- dimensional semiconductor systems. 
All the steps in this Section are independent of the single-electron basis used. H would 
be diagonal in the basis of the system's single-particle eigenstates, i.e., t^ lnl = e m \5 mini . 
In this paper we focus on the coherent response and we neglect coupling with phonons, 
which would result in additional, relevant dynamical variables and new contributions to the 
response function.-^ The SOS and QP pictures should be equivalent also when dephasing 
is included. In that case, however, the theory becomes more complicated. For the sake 
of simplicity and transparency we restrict the following analysis to the coherent response, 
where we do not include phonons explicitly. Dephasing effects, necessary for a realistic 
description, will be simply introduced by adding imaginary parts to excitonic frequencies. 
To introduce the exciton representation we define electron-hole operators:— 

R = A c Rt = J A\ 

±j m — u m2 u mi} ±J m — u mi u m2i 

where we have employed shorthand notation for pairs of indices: m = (mi, m 2 ). Using these 
operators we construct an effective Hamiltonian H (see App. IA~|) : 

H = J2 h mnBlB n + Urr^BlBlBkBi - £ [e + (t) fi* m Bl + E (t) Mm 5 m ) . (17) 

mn mnkl m 

The Hamiltonians H and are equivalent in the single and double excitations subspace, 
which is relevant for the response to third order in E.— This transformation from fermion 
to exciton variables is crucial for our approach, since it allows us to view the electronic 
degrees of freedom as a system of coupled oscillators. The parameters of the transformed 
Hamiltonian H are given by: 

h m n = ^m 1) n 1 ^m2i2 ^m\,n2^ m i n i W mirrl2nin2 , (18) 



mnkl 



j-W A A A _L f ( 2 ) A A A 4- 

L mi,ki u m2k2 u nih u n2l2 ' b m2to u m\k\ u n\l\ u n2h> 



fW A A A 4- AAA 

i ni ) h Um l k l Um 2k2 u n 2 l2 ~r L n 2 ,l2 " m l fc l u m 2 k2 u n 1 l 1 

t/C 1 ) A A -I- T/( 2 -* A A 

v min 1 kih u m 2 k2 u n2l2 T" v mzmkih ">1«1 "l'l 



+ 
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The commutation relations for the B operators can be obtained using the elementary fermion 
anticommutators: [c mi , CfcJ = 5 m i,fci- Within the subspace of |0) and B\ |0) states (i.e., the 
ground state and single excitations), we geiM 



B m ,Bl 



$mn 2 ^ l~ > mn pgBpBg , (1^) 



where 5 mn — S mini 5 m2n2 and 

Vmnpq &m\q\ ^p\n\ ^ni2P2 ^«292 ~^m2q2^P2fi2^m\p\^n\q\- (20) 

Eqs. (fT9l) and (I20I) are obtained in a similar way to (flTl) and (fl8|) . Terms with additional 
-Bj-Bj pairs (e.g. B^B^BB) are neglected in ( 1T91) . because they would introduce corrections 
higher than O (E 3 ) to the nonlinear response. Note the symmetry V mnpq = V mn qp- 

Using Eqs. (|T7|) and (fl9!) we obtain the nonlinear exciton equations (see Appendix [B]) for 
single-exciton variables ( I? m ) :-^>2£>^ 



d ( B m 

i X df 1 = h ™ (Bn) - »* m E + (t) + Vmnki (B n y (B k B^ (21) 

n nkl 



V 

npq 



+ 2E + (t) Vmnpq ( B n ) ( B q ) » 

where V is given by 



Vnmpq ^U nm pq 2 ^ nmlp^lq 2 ^ "Pnmkl Uklpq ■ (22) 

I k,l 



Here Y mn = ( B m B n ) are two-exciton variables. The Heisenberg equations give: 

E - ^ + W ((*») ^ + (^m) /<) (23) 



dY 



+ 2£+ (t)Y'Pmnkl(Bk) Vl, 

k,l 

Calculating the optical response by numerical integration of these equations^ 1 ^ is 
straightforward but numerically expensive. An alternative, more tractable approach, which 
further provides a better insight into the nature of the response, is to integrate the equations 
formally using one-exciton Green's functions G (t) and exciton scattering matrix T (t). The 
scattering matrix depends on quasiparticle statistics through the V matrix (Eqs. IB4| IB5j) as 

11 



well as on exciton-exciton coupling. This results in closed quasiparticle expressions for the 
3rd order contributions Si, Su and Sni to the response function (for details see Appendix 
□ and Ref. 3) 

si QP) (n,n,T 2 ,T 1 )= (24) 

T 43 r'J 

- 29 (r 43 ) 9 (r 32 ) 9 (r 21 ) Mn 4 M»n, / dr" dr' s x 



r^n'^'gni (j" 7"s) G n4n ^ (r^) G^ n3 (r 43 - r") G n ^ n2 (r 42 - r") G*/ ni (r 4i - r^) , 

where r 43 = r 4 - r 3 , etc. and G mn (t) = -%9 (t) [exp {-iht)\ mn . 

The response functions for the other phase- matching directions can be derived along the 
same lines. We get 

S < ff P \u,T 3 ,T 2 ,r l )= (25) 

-29(t 43 )9(t 32 )9(t 21 ) J2 Vn^VniHn, J K j K J2 * 

n 4 ...nx ^ n ' 4 ...n[ 

^n'^n^n'^s ~ r s)Gn 4 n' 4 (r' s )G n ' 3ns (r 43 - ?f )G*,„ 2 (t 42 - 7^)G n / m (r 4 i - r"), 

and: 

S?u\n,T Z ,T 2 ,n) = (26) 
-2£(t 43 )#(t 32 K(t 21 ) J dr'J J dr' s x 

n4-.n1 ^ r^...ni 

r„4„^ n ^„i(r" - ^)G , „ 4n ; i (r s ')G'^ ri3 (r 43 - r^)G n ^ 2 (r 42 - r")G n ' ini (r 4 i - r"). 

Just as in the SOS case, time translation symmetry implies that these response functions 
only depend on the three pulse delays t 3 ,t 2 ,ti. Eqs. (I24H26I) will be used next to connect 
the QP and the SOS pictures. 
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IV. CONNECTING THE SUM-OVER-STATES AND THE QUASIPARTICLE 
PICTURES 



We first recast Eqs. (I13H15I) using Green's functions (in all expressions t\ > 0, t 2 > 0, 
* 3 > 0): 

s (sos) = _ ^-Q^ h + h + h )fi-g{h)(i + G{t 2 )(i + ) (b3) (27) 

V G f (ti + t 2 )p,+&(t 3 )fi,-G(t 2 + t 3 )A + ) (c3) 
[jir& {h) /i + Gt(t 2 + t 3 )p,-G(t 3 )p,+) , (cl) 

Sf° s) = - (A"G(t3)A + G(t 2 )A"G ! (ti)/i + ) (a2) (28) 

V^te + t 3 )A~£(t 3 )A + G(ti + t 2 )A + ) (b2) 

A"G t (t 2 )A + G t (t 3 )A"G'(t 1 + t 2 + £ 3 )A + ) , (c2) 

S) = - (A~G(t 3 )A - ^(*2)A + G'(*i)A + ) (al) (29) 

V G^A'^te + t 3 )A + G(ti)A + ) , (bi) 



Here G(t) = —i6{t) exp(-iHt) and G^(t) = +i9(t) exp(iHt) represent the retarded and 
the advanced Green's function respectively; G, G and £y describe the evolution within the 
ground-state, single-exciton and double-exciton blocks of the Hamiltonian (Eq. ITT)) respec- 
tively. We also set the ground state energy e g to zero. 

Our goal is to show the equivalence of the QP and SOS pictures by deriving Eqs. (I24H26P 
from Eqs. ( T2TH29T) . To that end we adopt a harmonic reference system of noninteracting 
quasiparticle and expand the SOS response in anharmonicities. Harmonic oscillators are 
linear, and their nonlinear response vanishes identically^^ 2 .*^, as can be easily seen from 
the Heisenberg equations of motion. This means that the various Liouville space pathways 
for all nonlinear response function interfere destructively. Exploiting this property in the 
following derivation, we show that the quasiparticle physical picture has built-in cancellations 
in the reference harmonic system. 

We shall use the Dyson equation for the two particle Green's function, also known as the 
Bethe-Salpeter equation 33 

g (uj) = g H + g (co) r (cu) g h , (30) 
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or in the time domain: 



G(r) = Go(r) + [ T dr' f dr"Go (r - r') Y (r 1 - r") Q Q {r"). (31) 
Jo Jo 

Go is taken to be the Green's function of a doubly excited, harmonic system. It can be 
factorized into the product of a single-exciton Green's functions 

Go (''~)nin2n3n4 ^G nin3 {r^G n2nA (t) . 



The exciton scattering matrix V is defined by Eq. (J3TJ). 

Let us start with the Sj technique and show the equivalence of Eq. (|27|) to (1241) . The 
second and third terms of Eq. (J57j) (diagrams (cl) and (c3) in Fig. [H]) are purely harmonic, 
independent on the quasiparticle interactions. This is a direct result of the ordering of /Lt 1 * 1 , 
whereby the system only evolves in the ground and first excited state. Exciton-exciton 
interactions influence the evolution only in the second excited manifold. The first term 
in Eq. (l3Tj) . i.e. Go, represents harmonic evolution in the two-exciton manifold. Thus the 
first term of Eq. (1271) with G replaced by Go must cancel the other two terms, because 
the nonlinear response of a harmonic system vanishes. Substituting the second term from 
Eq. f[3"Tj) in Eq. (|2"TI) we obtain a single term for Sj: 

s (sos) = _ e 9 ^ f ia dr , r dr „ x (32) 

Jo Jo 

+ h + h)(i-Go(h - r')Y{r' - T")Go(r")fi + G(t 2 )^ 



The equivalence of the Eqs. ( |32|) and ( |24j) can be directly seen using the diagrams shown in 
Fig. ([H]). In these diagrams the scattering matrix V is represented by dashed regions. Note 
that G (t) G ] (t) = 6 (t) exp(-iht) exp(iht) = 6 (t). The QP diagram in Fig. © is obtained 
from the SOS one by changing the integration variables r' = t% — r' s and r" = — r'J. 
This completes the derivation of the QP expression for Si (Eq. |24|) starting from the SOS 
expression (Eq. |2T}) . 

Su can be calculated similarly. By combining Eqs. ( J28l) and (|3B the same type of 
cancellation of harmonic terms yields 

S (sos) = _ e ^ e ^ e ^ j tA dr , r dr „ x (33) 

Jo Jo 

'ft~G^ (h + t 3 ) (i-Go (h - t') T (r' - t") Go (r") ft + G(h + t 2 )tf 
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Eq. (J33]) is identical to Eq. ([25} as illustrated in Fig. (fTUj). 

We finally turn to Sin, (Eq. [29]). Using again the Bethe Salpeter equation (13T]) and the 
fact that terms that only depend on Q must cancel (harmonic reference), we get 

S { iu S) = -0(t 3 )0(t a )0(t 1 )x (34) 
jT dr' £ dr" (p,-G(t 3 )p,-g (t 2 - r')T(r' - r")^ (r")A + G(t 1 )A + ) 

+ f 2+< V f rfr' / (/i-Gt(t3) / i-4(t 2 + t3_ r ')r(r / -r / ')4(r // )A + G'(ti)^ 
Jo Jo x 

The equivalence of QP (Eq. [26]) and SOS (Eq. [29]) expressions can be shown as follows: the 
two terms in Eq. (1341 are labeled (SOSa) and (SOSb). The term (SOSb) can further be 
split into two terms (SOSbl) and (SOSb2), the first corresponding to r' < t 2 , the second 
to r' > t 2 (Fig. [H]). (SOSbl) is identical to (SOSa), but with opposite sign coming from 
&(t 3 ). Only the second term (SOSb2) remains, and it is equivalent to the (QP) diagram. 
We thus obtained Eq. fl26]) from Eq. ([3] 



V. 2D CORRELATION SIGNALS 

2D signals are displayed as correlation plots obtained by the double Fourier transforms 
of the various signals.— We shall denote the frequencies conjugate to the pulse delay times 
t\,t2 and £3 by fl±, fl 2 and f^. Starting with Eq. ([TT]) . and deleting some inessential factors, 
we obtain the induced polarization, which depends parametrically on the delay times t\, t 2 
and t 3 : 

P s {t 3 , t 2 , h) = S s (t 3} t 2 , ^y^ 1 "^ 2 "^ 3 " 3 ) V (Al " 1+A2 " a) V Alu;itl _ ^ 
Specifying the three possible signals by a proper choice of A factors we obtain: 

Pn{h,hM) = S//(t 3 ,*2^i)e i(wi - , * +ws)ta e < ( Ul - < *) ta e <Wltl , 
PmihMM) = S//i(t3,i2,ti)e ,(ui+urU3)l3 e !(ui+U2)l2 e" 111 . 

The 2DCS for Pj and Pn is defined as 

poo poo 

P a {n 3 ,t2M= dt 3 dhP a {t 3 , * 2 ,*i)exp {itt 3 t 3 + a = 1,11 (36) 

Jo Jo 
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For the SOS picture we use the expansions in eigenstates given by Eqs. (1Elj) . (1E2I) and (1E3I) . 
The QP expressions for Pj , Pjf P ^ and Pjfi are obtained along the lines presented 
in App. O Dephasing is introduced phenomenologically by adding a decay rate 7 to the 
Green's functions. We thus obtain Pj 

P} SOS) (n 3t t a ,n 1 ) = (37) 

% Y Vge^ge^ge^geJ^-^l + ^gih^gih^e^S ~ + ^2 + W 3 ) 
e2,ei 

+ i ^ HgetVleiVgex^teii-^l + ^l) J e 2 (*2)4i fo^ei (^3 ~ £*>1 + ^2 + W 3 ) 
e2,ei 

- « XI ^/Mei/MgeiMsea^C-^l + ^l)4* 2 C^)^ (^)^>e 2 (0 3 - + W 2 + W 3 ), 
e2,ei/ 

P / (QP) (fi 3 ,t2,^i)= (38) 

- 2 X fl ei fl* e3 fl* e2 fl ei r ei {t2)Ie 2 {t2)I* ei {-^l ~ Wi)/e*(fl3 ~ ^1 + ^2 + W 3 ) 

ei,e2,e3,e4 

X r e4ei e 3e2 (fi 3 - + W 2 + ^3 + £ej + «7ei)£oe 3 e 2 (^3 ~ ^1 + ^2 + ^3 + £ei + HeJ- 

The Green's function Fourier transform is defined as G(w) = / dtexp(iut)G(t) [and G(t) = 
/ f*exp(-iw*)G(w)]. We have 

J e (cu) = ^e|G(w)|e^ = (cu-5 e + z 7e )- 1 , (39) 



^oe 2ei M = (e 1 e 2 



yoM 



eie 2 



1 



(40) 



io - e e2 - e ei + z (7 e2 + 7 ei ) 
Eq. (|40l) is obtained by transforming Qokijr (t) to the single-exciton basis and performing the 
Fourier transform. We also define 

Fab{t) = -i0{t) exp (z (e 6 - e ) i - + 7 6 ) t) , 

Fab(uj) = (u - e a + e b + i7„ + i^) -1 . 

Similarly we obtain for P77: 

p r ( f os) (n 3 ,*2,fii)= (4i) 

X Vl^geiVgextfgeJexiSll + ^l)I* g {t 2 ) I g {t 2 ) I e , 2 (Q 3 + W X - U; 2 + Wg) 



77 



e2,ei 



X MseiMge 2 Mse 2 A*g e J ei (fil + UJi)I* 2 {t 2 )I ei {t 2 )I ei {Q 3 + U>1 ~ ^2 + W 3 ) 



e2,ei 



+ * ^/^/^e, Mg e / ei (^l + CJi)/ e * (t 2 )/ ei (t 2 )^7e 2 (^3 + ^1 - ^2 + W 3 ), 

e2,ei,/ 
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p!f P) (Q 3 ,t 2 ,Q l )= (42) 

- 2 ^ ^e i K i ^e 2 ^lJl 2 {h)Ie 1 { t 2)h 1 {^l + Wl)/ e4 (^3 + ^1 ~ ^2 + W 3 ) 

e4..ei 

X r e4e2e3ei (fi 3 + - LU 2 + ^3 + £ e2 + He 2 )Goe 3 ex(^3 + ^1 ~ ^2 + ^3 + ^e 2 + «7e 2 )- 

The P/// 2DCS signal is denned as 

poo POO 

Pni{n 3 ,^2,ti) = / dt 3 / dt 2 Pin(t 3 ,t 2 ,t 1 )exp{in 3 t 3 + in 2 t 2 }. (43) 
Jo Jo 

This yields: 

pg os \n 3 ,n 2 ,t 1 )= (44) 

-i fj, gei iJ, ei fiJ,* e2f iJ,* ge2 I e2 (t 1 )If(£L2 +wi +a; 2 )4i(^3 + ^2 -^3) 

e2,ei,/ 

+ i ^ ju flei /i ei/ /x* 2/ //* e2 J e2 (ti)/ / (n 2 + a;i + a;2)^/e 1 (fi3+^i+^2-^3), 

e2,ei,/ 

P / ( / ( ? P) (fi 3 ,fi 2 ,t 1 )= (45) 

- 2 2J ^MeaMeaMe^ei^l)^!^ + ^1 + ^2 ~ W 3 )/* 3 (f2 2 ~ ^3 + ^i) X 

e 4 ...ei 

[r e4e3e2ei (fi2 + ^1 + ^2)^0e 2 ei(^2 + ^1 + ^2) 

-r e4e3e2ei (fi 3 + + W 2 - ^3 + ^3 + *7e 3 )£oe 2ei (^3 + + ^2 ~ ^3 + £ e3 + ij ea )] • 

Both p( 5 ° 5 ) and p( QP ) depend on the single-exciton energies. However, the SOS expres- 
sions contain two-exciton eigenenergies (Ef) explicitly, while the QP counterparts contain the 
scattering matrix T instead. The equivalence of the two representations has been established 
in Sec. [TV] Eqs. ( 13711451) constitute our final expressions for the various 2DCS signals. In 
this form they may be readily used in numerical simulations. The SOS expressions (Eqs. [371, 
|4~T1 and l4"4l) were recently used to survey the various possible resonances and cross-peaks in 
2DCS of semiconductors.— 



VI. DISCUSSION 

The quasiparticle representation is obtained using the Heisenberg equations for the ex- 
citon oscillator variables. These equations form an infinite hierarchy involving successively 
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higher numbers of excitons.-^i The hierarchy may be truncated, depending on the observable 
of interest. For instance, the absorption originates from single-exciton creation/annihilation. 
Only single-exciton variables should then be considered, and exciton-exciton interaction 
terms may be neglected. The nonlinear response depends on the exciton interactions, thus 
single- and double-exciton variables need to be treated explicitly. The two coupled NEE 
equations (|2lT 123]) describe the third order response. These equations are exact in the ab- 
sence of dephasing. When dephasing is included by adding linear coupling to a phonon bath, 
two additional variables ( B^B) and ( B^BB) must be included in the NEE to describe the 



cle. For some techniques the present equations provide a good approximation even in the 



The quasiparticle approach avoids the explicit calculation of multiple exciton states: their 
influence is represented by the scattering matrix, which can be calculated provided the exci- 
ton interactions are known. We have shown how the quasiparticle expressions for the various 
third order techniques, ordinarily derived by solving equations of motion, can be obtained 
directly from the sum-over-states expressions by employing the Bethe-Salpether equation. 
These expressions explicitly contain the two-exciton Green's functions and have many inter- 
fering terms with large cancellations,— which complicate their numerical implementation. 
In the QP picture, on the other hand, these interference effects are built-in, considerably 
simplifying the expressions for the nonlinear response. 35 

The interpretation of 2DCS signals using the SOS expressions is straightforward.— In 
the k/ technique one-exciton coherences are observed during t\, and the coherences between 
excitons and biexcitons are observed during £ 3 . Thus the 2DCS shows peaks along Qi and 
f2 3 corresponding to these resonances, km shows biexciton resonances along the Vl 2 axis, 
this technique is known in NMR as double-quantum coherence. We have established the 
connection between the SOS and the QP pictures by using time-ordering on the Keldysh- 
Schwinger loop, which only maintains partial time ordering in real (physical) time. 

The Hamiltonian H (Eq. ITT)) can describe several microscopic models other than the Wan- 
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nier excitons considered here (Ht)- Vibrational excitations (soft-core bosons) and Frenkel 
excitons (hard-core bosons, Paulions) in molecules can be mapped into the same model.— isjSj&I 
The equations of motion for these other systems are similar, but not identical, because of the 
different commutation relations (QP statistics). Eq. ffl9l) provides a unified description for 
all of these systems, by specifying the proper commutation rules:- 1 ^ for bosons V mnpq = 
and for Paulions Vmnpq = 8 mn 5 mq 5 np . These expressions for V may be substituted into our 
final expressions for the response functions, where they only affect the exciton scattering 
matrix, which in the frequency domain reads (App. [D]) 

r(u) = (i — vg {u)y l vg Q (u) (i - v) g \uj) - vg \uj), 

where V is given in Eq. (l22l) . Q is the free two-exciton Green's function (App. [Bj and I 
is the tetradic identity matrix. The nonlinearity of the system depends on QP interactions 
as well as non-boson statistics; both enter through T. For noninteracting bosons, where 
U = V = 0, T vanishes and so does the nonlinear response. In Appendix [F] we present V for 
bosons and Paulions. 

We have used the symmetry P mnpq = P mn qp in our derivation. Since the boson commuta- 
tion relations are simpler than for Fermi or Pauli operators, a considerable effort has been de- 
voted to mapping the original problem with complicated commutation relations into a boson 
picture.— 1 ^ The resulting boson Hamiltonian contains additional interactions which com- 
pensate for the statistics. For instance, the Frenkel exciton Hamiltonian for Paulions may be 
mapped into an anharmonic Hamiltonian of bosons with quartic couplings. Bosonization^ is 
very convenient for describing exciton scattering: the response functions derived for bosons 
can be applied for arbitrary operators, provided we modify the Hamiltonian and express it 
in terms of boson operators. 
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APPENDIX A: EXCITON REPRESENTATION OF THE TWO-BAND HAMIL- 
TONIAN FOR FERMIONS 

By construction, the Hamiltonians H (Eq. (1171) ) and H? (Eq. (1161) ) are equivalent only in 
the physically relevant space of single and double excitations. This is sufficient to calculate 
the response to third order in the field E(t). H may be constructed using the following 
rules: 

• since the Hamiltonian (fl6l) conserves the number of excitons, it should only contain 
products with equal number of B^ and B operators (except for the Hi term, which 
does change the number of excitons) 

• a term B^B^ . . . B\ B^B^ . . . Bb p gives zero when acting on states with less than p 
excitations and only affects manifolds with p excitations and higher. 

The parameters of H can be obtained as follows. First we note that no con- 
stant term k should be added to (jT7|) . since it would yield: (0 \k\ 0) ^ 0, while 
(0\H T \0) = 0. The Y. m i,n2,k2A Wm ^ l ^ c ^i d U d k2 c h term of h t can be written directly 
as ^ mi , m2 ,n 1 ,n2 W ^irn2n 1 n 2 cl ni dl n2 d n 2C nl = Y, m ,n W mnB ] m B n . Also the term describing the 
interaction with light can be obtained directly. Using the second rule given above we imme- 
diately see that no terms higher than B^B^B^B^ are necessary in the sub-space defined 
by functions |0), Bj |0) and B}Bj |0). We thus obtain the form given in (fTTl) . We next 
calculate, in this sub-space, matrix elements of H, and compare to matrix elements of Ht- 
In this way a one-to-one correspondence of the parameters of H and H? can be established. 

Additional terms must be included in H in order to describe higher order response func- 
tions. This can be done using the same rules. 



APPENDIX B: THE NONLINEAR EXCITON EQUATIONS 

The Heisenberg equation of motion (NEE) for the Hamiltonian ( TT7|) reads: 
d(B n 



dt 



/ J h nm (B m j — H* n E + + 2j Vnmpq (^BpBg) (Bl) 



rnpq 



rnpq 
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Here we invoked RWA and used the notation of Eq. ([I]). Employing (1201) we see that V nmpq = 
Vnmqp, so the last two terms in Eq. flBTj) can be recast as: 2E + ^2 mpg V n mpq ^4r4) fj,*. 
We now make the following factorization: 



4n4 



B m ) (B p ) and (4,44) = (4.) (44/ 



(B2) 



which is exact for pure states when dephasing is neglected^ and is a good approximation 
in the absence of incoherent exciton transport. Eq. (1B1I) then yields the Eqs. (l2Tj) and (1231) . 
where 



h 



(Y) 



%,kl ~ fimkhnl + bnlhmk + Vmnkl = h + V. 



(B3) 



(i) 



We next expand the EOMs in orders of E. Using B™' for ( B m } we obtain: 
^/w4 1} - »* m E + (t), 



.dBm 
dt 



dY {2) 

. w i ran 
I : 

dt 

.dB^m 
1 dt 



E h %U Y * - E+ (*) + + 2£ + (t) E r mnkl B^^t 



kl 



k,l 



The Green's function (tetradic matrix) for F^ 2 ^ is £ (t) mnkl = 
thus 



-i0(t) [exp (-z/^t)] 



mnkl ' 



£ £w* (t - r) (r) (4 1} (r) M * + B« (r) tf) - 2 £ T^flW (r) /** 

/CxD /"CO 
dr / dr' £ </ mnjW (t - r) x 
CXD J —CO ,„1„ 



(r - r ) /4 + G fea (r - r) /x* - 2 £ 



,,„ v r - r) ^* 



(r) (t') 



We also define the zero-order tetradic Green's function Qq for for the case V = 0, i.e. 
Gomnki (t) = —id (t) [exp (— iht)\ ,,, it will be used later. Bm is given as 



/oo 
G mn (t' - t) x 
-oo 



£ W#>* (t) y (t) + 2£+ (t) £ K fe p ? 4 1} * (*) 4 1} (*) 



dt. 
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This expression can be simplified using the symmetry Gkifg = Gkigf- At this point we 
introduce the tetradic exciton scattering matrix T defined as: 

T (u) Go (w) = VG (u) (I — V)— V, (B4) 

which in time domain can be written as (see App. ID]) : 

oo 

VG (t - t) (I - V) = V8 (t - t) + J dr x Y (t - t x ) G (n - r) , (B5) 

— oo 

where the tetradic identity matrix If g j r = 8fj$gr- Since G (t — r) ~ 9 (t — r) is retarded, T 
must be retarded as well, i.e., r (t — Ti) ~ 9 (t — T\). To proceed further we take advantage 
of the factorization: 

Gokijr(t)=iG kj (t)G lr {t), (B6) 

which can be easily shown in the single-exciton eigenbasis. After a rearrangement of terms 
we obtain: 




B$ (r 4 ) = -2 / / / / / dT'dT 2 dT 1 dT"dT 3 



G «4< ( r 4 - r ") r^n;^ (^" - G;- n3 (r" - r 3 ) x 

G n>2 (r' - r 2 ) G n[ni (r' - n ) 6 (r 2 - n) <XiMn 3 £ + (r 2 ) E+ (n) E~ (r 3 ) 
The 3rd order polarization is 

(T4 ) = y: + « )f ) 

114 

= -2 J J J dr z dr 2 dT X ^ Mn 4 Mn3^ 2 M^ 

124,113,112,111 

9 (r 2 - r x ) J dr" J dr' ^ T (r" - r') x 



n 1 ,n2,rtg,n 4 



G„ 4 „^ (r 4 - t") G* n ^ n3 {t" - r 3 ) G n > 2n . 2 (r' - r 2 ) G n > ini (r' - n) X 
£T (r 3 ) (r 2 ) (n) + complex conjugate, 



where we used r n > „/„/„> (t) = (t). 
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The above expression is finite only for t 2 > T\. Hence there are 3 possible intervals for r 3 , 
that define three contributions to the third order response function S^ p ^ 



00 00 



(r 4 )= J dr 2 j dnO^-T^x (B7) 

—00 —00 

T\ T2 +OO \ 

j Sf P) dr, + Jsfi p) dT 3 + J Sf I PdTAE-{n)E+{T 2 )E+{T 1 ) + c.c. 

-00 n t 2 / 

This definition of Sf P \ 5^ P) and S$P is used in Appendix (0 to obtain the Eqs. (J24H26D- 

APPENDIX C: RESPONSE FUNCTIONS OF QUASIPARTICLES 

For calculating each of the contributions to Eq. (1B7I) we need to switch to a different set 
of time-ordered variables. For S\ QP ^ we set: T\ — > r 2 , t 2 — > r 3 , r 3 — > n: 

-+00 



/+00 r+00 
dr" dr' > 

r n in^»»i ( T " ~ T ') G "4< ( r 4 - 7"") G*- (r" - T X ) G n ^ n2 (r - T 3 ) G n ' ini (r - T 2 ) 



where r 4 > r 3 > r 2 > T\. 

Substituting t' s = t 4 — r" and r" = r 4 — r' and exchanging the dummy indices ni — >■ 
n 2 , n 2 — > n 3 , n 3 — ► ni (same for primed indices) we obtain Eq. fl24|) . Integration limits for 7^ 
have been limited by G nAn < (t' s ) and T (r" — t' s ), while for r" by G n ^„ 3 (r 43 — r"). Eqs. (|25|) 
and f[2T)j) are obtained in a similar way. 

Eq. ([MP can be simplified considerably by performing the double time-integrations ana- 
lytically. We first express the exciton Green's function G (r) and T (r) in the one-exciton 
basis ip e , defined by: 

^ ] hmn^en = ^e^emi (CI) 
n 

where /i mn is given by Eq. (fI5|h The energies e e define the lowest optically-excited manifold 
of the system, i.e., single excitons. In this basis we express the time and frequency-domain 
one-exciton Green's functions: 

G mn (r) = "Pernle (r) h (r) = -id (r) exp (-ze e *) , (C2) 
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where we introduce dephasing via e — > e — ry e . The tetradic exciton scattering matrix is 
given by 

r m4 r)i3jn2r?ii ( r ) — ^ ^ V'e4m4 "^637713 r e4e3e2ei (t) V^m^eimi • 
el...e4 

and the transformed dipole matrix elements /x e9 are fi e = Ylim Vmipem- Using these quantities 
we express the S ( j QP) in the single-exciton basis: 



S\ QP \r 4 , r 3 , r 2 , n) = -20 (r 43 ) (r 32 ) (r 21 ) ^ (C3) 



ei...e4 

T43 

s-l e4e 3 e 2 ei l 



'0 ./0 

x ^ 4 (T s Ve 2 (T43 - r' s ')I ei (r 42 - r' s ')I* 3 (r 41 - t' s ). 

We next introduce T (t) = J ^e~ iuJt T (u). Since the response function depends only on the 
pulse delays and not on the absolute times, we denote S , /(r 4 , r 3 , r 2 , T\) = Si(t 3 , t 2 , ti), where 
^3 = t 43j ^2 = r 32 and = r 21 . We perform a Fourier transform with respect to the first and 
last arguments. We thus obtain 

5f p) (fi 3 , t 2 , no = -2*6/ (t 2 ) ^ Me 4 M: 3 M: 2 Me 1 /; 1 (-no n i ea n h A m x m 

el...e4 

y dul^^ (w) e0e 3 e 2 (w) ^ (w ~ ^s) , 

The uj integration can be performed by noting that 

which is obtained from Eq. ( 1B4I) by noting that (a;) has poles only at two-exciton energies, 
and that 2e — 2i r ) is a good approximation for two-exciton energy and dephasing rate. Hence, 
if we close the Cauchy integration path in the positive half-plane, there will be only a single 
pole at uj = Q 3 + e ei + ij ei as seen from fl39l) . This finally gives 

Si QP) (fi 3 ,t 2 ,fiO = -2 ^e 4 M: 3 M: 2 Me 1 /: i (*2)/e 2 (t2)/: i (-nO/e 4 (n3) (C5) 

e 4 ...ei 

X r e4eie3e2 (fi 3 + £ ei + Hl)G e 3 e 2 {^3 + £ei + *7l)> 

To account for carrier frequencies tui, cj 2 and u; 3 appearing in the polarization (Eq. [35]) at 
this stage, we can perform the substitution — > Qi + oj\ and Q 3 — > —Qi + uj\ + u> 2 + u> 3 . 
In this way we obtain Eq. fl38l) . Eqs. fj42l) and (j4"5]) are derived similarly. 
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APPENDIX D: THE EXCITON SCATTERING-MATRIX 

In order to use equations (12411261) for calculating the quasiparticle response function, we 
should calculate the scattering matrix T. We first write Q (u) and Go (lo) in an operator 
form 

Go (w) 



lu — h 



oj - hPO oj-h-V 
where h is defined in (IB3j) . The Dyson equation then reads 

Q = Go + GoVG (w) = Go + GoVGo + GoVGoVGo + • • • , 
which can be recast in the form 

vg = (i - vGoT 1 vGo- 

Using Eq. flB4j) . we obtain: 

TGo = (I — ^So)" 1 VGo (1-V) - V, 
which results in the final expression for T 

r = (I - VGoY 1 VGo (i - V) Go 1 - VGo 1 (Dl) 
The l.h.s. of Eq. (IB4j) can be expressed as a convolution: 

y dr' y dhT (h) g (r - t{) exp (W) . 
The r.h.s. can be written as 

V J dr'Q (r) (I - V) exp (W) - J dr'VS (r') exp (W) , 
note that V is independent on r' or cj. Since l.h.s. =r.h.s. for any ui, we must have: 

y d^r (t x ) £ (r' - tO = ^ (r') (I — P) — W (r') . 
Substituting r' — > t — r and ti — > i — n we obtain Eq. (1B5I) . 
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APPENDIX E: SOS EXPRESSIONS FOR THIRD ORDER TECHNIQUES. 

Upon expansion in the eigenstates for the exciton level scheme shown in Fig. [3] we get 

e,e' 

(fl-fl-fl + fl + ) = J2J2(V g e>Ve>f»jeVtg) ■ 

e,e> f 

Expanding Eqs. ([13] - USD in the eigenstates, we obtain the sum-over-states expressions 
for the third-order response functions: 

Sf° 5) (t 3 ,t 2 ,tl) =i i e{h)d{t 2 )d{t 1 )Y,» 9 e>»e>gVgeVegi: i {tl)Ie{h) (El) 

e,e' 

Hge' He'gHgeHeg 

e,e' 

e,e> f 
e,e' 

H'ge' He'gHgeHeg /*, (t 2 ) / B (tj + t 2 + t 3 ) 

e,e' 

Mge'Me'/M/eMeg-^e' (*2 + ^) ^/ (* 3 ) (*1 + ^ 2 ) , 

e,e' / 

Sjff^Ml) = ^9{h)9{t 2 )e{t X ) J^J^^e'Ve'fflfeflegle' (* 3 ) If (t 2 ) J e (t X ) (E3) 

e,e' / 

-< 8 e(*3)^ a )f(«i)EE 

V>ge> He' fl*f elegit' (^3) -fy (*2 + ^3) (^l) , 

e,e' / 

where I e (t), defined in Eq. (]C2p . is the Green's function in the single-exciton eigenstate 
basis. Eqs. flBTHTOD immediately follow by substituting Eqs. (IE1IIE2IIE3D in Eq. ([36]) and 

APPENDIX F: QUASIPARTICLE PICTURE FOR SOFT-CORE AND HARD- 
CORE BOSONS 

In this Appendix we apply our QP expressions to two other types of quasiparticles with 
different statistics. These two examples demonstrate the generality of our approach discussed 
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briefly in Sec. IVll 

We first consider the Hamiltonian of a system of coupled anharmonic oscillators (soft-core 
bosons) : 

mn mnkl 

where B^ and B n are boson creation and annihilation operators with commutation 
[B m , B£] = 5 mn h mm is the fundamental transition energy of the mth oscillator, while h mn 
is the coupling between the mth and the nth oscillators. U mn ki is the anharmonic coupling. 
This Hamiltonian has been used to describe infrared nonlinear spectra of proteins.— ^ 

For this model the scattering matrix can be obtained from (1D1I) by putting V = 0. In the 
site representation it reads: 

r = (i-vg y 1 v, 

here T is a tetradic matrix, V = 2U and Go (^) is defined in Eq. ( 1301) . 

We next turn to electronic excitations in molecular aggregates or crystals with weakly 
interacting molecules. These are described using the Frenkel Exciton Hamiltonian. If the 
excited-state absorption frequency of each molecule is well separated from the ground state 
absorption, the excitations can be modelled as coupled two-level systems.— 1 ^ The Hamilto- 
nian is 

H h mn B^ m B n . 

mn 

The nonlinearities are now hidden in the statistics of exciton creation (J3+) and annihilation 
(B n ) operators. These are bosonic for different oscillators (units) and fermionic for the 
same oscillator. Their Pauli commutation relation is [B m , = 6 mn ^1 — 25^5 n j. The 
commutation relation ensures that two excitations are not allowed to reside on the same site 
(hard-core bosons). The scattering matrix in this case is given by: 

I" ^ mnkl $mn$kl^mn j 

f = -Q(io)- 1 , 

and 

Gmn (^0 &mmi$nn\Gomm\nn\ y^J ■ 

This form of the exciton scattering matrix was recently successfully applied to study molec- 
ular chirality induced signals in molecules.— It can be obtained from Eq. fIDlj) in the limit 
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U = 0. All QP-statistics effects (Paulion commutation relations) are included in the relation 

V nm pq 2 T~ > nmlphlq (Ecj- 122 j) and l~ > nmlp ^nm^nl^mp- 
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FIG. 1: The sequence of light pulses in a time-domain Four Wave Mixing Experiment: the pulses 
are centered at times t±, T2, T3, while the delays are t±, t<i and t%. (The latter are sometimes 
denoted as r, T and t.) The signal is generated in the k$ direction. 
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FIG. 2: Diagrams representing the four partially time-ordered terms contributing to the third-order 
polarization (Eq. [5]). 




FIG. 3: Energy levels of the exciton model. 
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FIG. 4: Loop diagrams representing the partially time-ordered terms contributing to the third- 
order polarization [Eq. ([8])] within the rotating wave approximation. Arrows pointing to the right 
represent n + and arrows pointing to the left /x - . The diagrams are obtained by adding arrows to 
the interactions in Fig. ([2]). 
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FIG. 5: Diagrams representing the fully time-ordered terms contributing to the third-order polar- 
ization within the rotating wave approximation (Eq. (J9|)). Tj repesent the interaction times with 
the various fields. Arrows pointing to the right (left) represent /x + (/x~). Time variables in loop 
diagrams (b) and (c) on the left are ordered in the loop. The other open diagrams are fully ordered 
in physical time. 
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FIG. 6: Feynman diagrams for the kj technique (Eq. |13|) . 
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FIG. 7: Feynman diagrams for the kjj technique (Eq. [T 
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FIG. 8: Feynman diagrams for the km technique (Eq. I15|) . 
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FIG. 9: Loop diagrams showing the equivalence of the Si expressions in the QP (Eq. I24D and SOS 
(Eq. 1321) pictures. Dotted, single and double lines show ground, single exciton and double exciton 
states' evolution respectively. Dashed region represents scattering matrix. 
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FIG. 10: Loop diagrams showing the order of time variables in the QP (Eq. [25]) and SOS (Eq. [33]) 
expressions for Sn- 
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FIG. 11: Loop diagrams showing the order of time variables in the QP (Eq. 126ft and SOS (Eq. 134ft 
expressions for Sin- (SOSbl) cancels with the (SOSa) term in (Eq. 134ft (not shown). The remaining 
diagram (SOSb2) is identical to the (QP) diagram with a simple change of time variables. 
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